Electrostatic interactions between discrete helices of charge 
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We analytically examine the pair interaction for parallel, discrete helices of charge. Symmetry 
arguments allow for the energy to be decomposed into a sum of terms, each of which has an intuitive 
geometric interpretation. Truncated Fourier expansions for these terms allow for accurate modeling 
of both the axial and azimuthal terms in the interaction energy and these expressions are shown to be 
insensitive to the form of the interaction. The energy is evaluated numerically through application 
of an Ewald-like summation technique for the particular case of unscreened Coulomb interactions 
between the charges of the two helices. The mode structures and electrostatic energies of flexible 
helices are also studied. Consequences of the resulting energy expressions are considered for both 
F-actin and A-DNA aggregates. 

PACS numbers: 87.15.-v 



I. INTRODUCTION 



Many important biological polymers are both acidic 
and helical in structure. Specific examples include DNA 
and F-actin, a key component of the cellular cytoskele- 
ton. When placed in water, each of the acidic subunits 
of these polymers becomes negatively charged. This re- 
sults in a discrete electrostatic charge distribution which 
follows the shape of the polymer's helical backbone. Be- 
cause these molecules each carry a large net negative 
charge, they will strongly repel one another under typ- 
ical conditions. However, through the introduction of 
multivalent counterions, cross-linking agents, or osmotic 
stress, the molecules can be condensed to a high den- 
sity 0, Q • This ability to induce aggregation of like- 
charged molecules is important for biological systems. 
Compacted DNA is found in cell nuclei, bacteria, and 
virus capsids, for example, while F-actin bundles play an 
important role in the processes of cell motion and division 
0,11, [1,0, Hi- At high density, the interaction potential of 
two neighboring molecules will depend strongly on their 
relative positions and on the many parameters specify- 
ing their helical structure. This sensitivity can result in 
various conformational phase transitions dependent upon 
lattice and individual helix symmetry couplings [1| . Un- 
derstanding the physical mechanisms behind these tran- 
sitions and how they are used for biological control has 
been one of the major goals of the literature considering 
this topic. 

To model these systems many previous theoretical 
studies have considered the interactions between contin- 
uous helices of charge. Such considerations led to in- 
teresting results which appear to rationalize various ex- 
perimentally observed phenomena. These include the 
possibility of a B to A conformational phase transition 
for condensed DNA bundles [lO, UM as well as more 
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general frustration-induced packing geometry transitions 

Although physical systems always consist of discrete 
charge distributions, theoretical studies of the interac- 
tion between discrete helices of charge have been com- 
paratively limited. Various studies have considered 
particular integer number of charges per turn systems 
[lOl . [3, [H, liSi but to the authors' knowledge, only one 
previous study has looked at the general discrete case 
[l7| . In that reference it was demonstrated that the in- 
teraction energy of two parallel helices is a discontinuous 
function in the angle between charges on a single he- 
lix. This remarkable result followed from the fact that 
helices with a rational number of charges per backbone 
turn, hereafter referred to as rational helices, may in- 
teract through modes which vanish for a pair of corre- 
sponding irrational helices. Here, the modes considered 
are terms in the series representation obtained in Ref. [ll| 
for the interaction between two general cylindrical charge 
distributions, explicitly taking into account the effects of 
both adsorbed and free counterions. Modes unique to 
rational helices were found to allow for energy reduc- 
tion. This discrete effect could possibly play a significant 
role in determining the conformational twist of interact- 
ing helices. For example, through the consideration of 
an idealized, perfect helical model, it was shown that 
this mechanism could provide a sufficient amount of en- 
ergy to allow for the twisting of B-DNA from its isolated 
value of 10.4 charges per turn to the 10.0 charges per turn 
observed in aggregate Although the same authors 

later argued that sequence dependent variations in twist 
should wash out this effect for B-DNA [1|, their earlier 
work indicates that the consequences of discreteness can 
be significant. Thus, in principle, the phase diagrams of 
these systems should depend not only on the smoothed 
out helical shape of their constituent molecules, but also 
on the symmetries of their discrete charge distributions. 

In this paper, we revisit the interaction between a pair 
of helical discrete charge distributions. A simple model 
is considered in an effort to focus directly on the sym- 
metries of the interaction. These symmetries allow us 
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to demonstrate useful relationships between the interac- 
tion energies of each of the three most commonly studied 
models for these systems. These include the discrete he- 
lical, continuous helical, and cylindrical models, the lat- 
ter being a model in which each molecule is replaced by 
a continuous cylindrical surface charge. This approach 
allows for a convenient decomposition of the discrete in- 
teraction energy which takes the cylindrical interaction 
energy as a base energy and then adds on independent 
correction terms associated with different aspects of the 
helical shape of the molecules. The rational discrete cor- 
rection terms, equivalent to the extra modes discussed in 
Ref. [ll], are considered in some depth. 

Algebraic analysis allows us to determine information 
regarding the phases and amplitudes of the Fourier series 
expansions for each of the correction terms. For exam- 
ple, an extra symmetry of the system is considered and 
is shown to cause a reduction in the energy benefit of ra- 
tionality for certain high-symmetry orientations. In ad- 
dition, each of the series are shown to be exponentially 
convergent. Truncation of these scries therefore allows 
for simplified, approximate expressions for the energy to 
be obtained which capture the basic structure of the en- 
ergy landscape observed numerically. Although we focus 
on the electrostatic interaction, it is later shown that the 
same expressions also apply quite generally and can be 
used to model other forms of interaction. We apply the 
effective energy expressions to the examples of F-actin 
and A-DNA to demonstrate how they may be used to 
make strong statements regarding the azimuthal inter- 
actions and frustrations in these systems. As discussed 
previously [l0| . reduction of azimuthal frustrations may 
help to explain why particular lattice packing structures 
and twist angles are observed in experiments. 

In addition to the static helix-helix interaction energy, 
we also consider the mode structures of both an isolated 
and an interacting pair of flexible helices. An instability 
is indicated by the single helix mode analysis for certain 
choices of the angle between charges. To understand this, 
we briefly consider the electrostatic energy of a single he- 
lix of charge. For a pair of interacting helices, we show 
that a gap is expected at long wavelengths between the 
oscillation frequencies of a rational system and those of 
a corresponding irrational system. Appropriate neutron 
scattering experiments may thus provide information re- 
garding the rationality of twist angles in bundled sets of 
helices. In principle, such experiments could directly test 
whether sequence dependent twist variations indeed wash 
out the effects of rationality in B-DNA systems. This was 
argued to be the case by the authors of Ref. ^ on the 
basis of angle variations inferred from x-ray diffraction 
data. Angle variation need not only result from base- 
pair sequence variation, however. Local twisting should 
also occur within these systems, resulting in an increase 
in the energy beneflt associated with discreteness. This 
suggests that further study of disordered, flexible sys- 
tems is required in order to rule out discreteness as the 
cause of twisting in B-DNA aggregates and elsewhere. 
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FIG. 1: The fundamental configuration considered. 



The present study provides a first step in this direction. 

The paper proceeds as follows: the mathematical for- 
mulation for the problem considered is presented in the 
following section; in section III geometric symmetry ar- 
guments are presented which provide the basis for the 
energy decomposition; section IV addresses the effects of 
symmetry on the phases and amplitudes of the Fourier 
components of the terms in the decomposition; section 
V contains the energy and mode analysis for azimuthally 
flexible systems; section VI explores the possible conse- 
quences of our simplified energy expressions for physi- 
cal F-actin and A-DNA aggregates; section VII contains 
concluding remarks; finally, an outline of the Ewald sum- 
mation technique applied to obtain quick numerical eval- 
uations of the energy appears in an appendix. 



II. MATHEMATICAL FORMULATION AND 
NUMERICAL EVALUATION 

The fundamental configuration of charges considered 
is depicted in Fig. [TJ Two identical, infinite helices of 
discrete charge lie parallel to one another, each with ra- 
dius r and separated by a distance d. The interaction 
energy of the two helices may be formally expressed as 
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E 



exp[ 



R 



(1) 



where the sum is over all pairs of charge, one taken 
from each of the two helices. Here we have assumed a 
Yukawa type individual charge potential, consistent with 
the Debye-Hiickel screening approximation , is the 
screening parameter, and the distance R between two 
charges is given by 



[d-r cos(0i - 0) -I- r cos(6'2 - 0)]^ 

+ [rsin(0i - 0) - r sin(02 - <P)f + a'[^i - ^2 + C]' 



2rd[cos(6'i -</))- cos(6'2 - (j))] 
2r 



+2r"[l - cos(ei -62)]+ a'[ei - 9^ + ( 



(2) 



3 



FIG. 2: The helical charge distributions, shown end-on. 



In the above, the parameters 9i and 62 specify the az- 
imuthal angles of the two charges being summed over, cj) 
is defined in Fig. [21 a is related to the helical pitch, and 
C describes a vertical or axial shift between the two he- 
lices. To begin, we assume perfect helices and the angles 
9i and 62 are written as 



h = ni4' 

h = n2tjj + At/j, 



(3) 
(4) 



where the rii are integers to be summed over, ^ is the 
azimuthal angle between adjacent charges on a single he- 
lix, and Alp is the azimuthal shift of the charges on the 
second helix relative to those on the first. For rational ip 
we write 



(5) 



with and relatively prime. Equation ([5]) indicates 
that each helix has charges for each m.0 turns of the 
backbone. It follows that the discrete energy is, in gen- 
eral, periodic in C with period 2Trm^,, in Aip with period 
i/), and in with period 2Tr/n^. 

All numerical values of the energy presented here are 
for the case of Coulomb interactions with no screening. 
That is, Qs was set to zero. This was done to simplify 
the expressions being summed over; while this somewhat 
removes the model from the physical systems considered, 
the numerically calculated energies retain the symme- 
tries of these systems. It is these symmetries which are 
the focus of study in this paper, and fortunately, these 
are quite insensitive to the individual potential functions 
used. The following two sections demonstrate how these 
symmetries allow for a characterization of the interaction 
energy. Generalization of these results will be discussed 
in section VI. 



III. GEOMETRIC SYMMETRY ARGUMENTS 

Consider the potential felt by a charge on one irrational 
helix due to the charges on a second, identical helix. Lo- 
cally the second irrational helix will look like a helix with 
a nearby rational number of charges per turn. Since the 
potential of the second helix is a continuous function in 
tjj, it follows that the potential energy of the first charge 



is the same as it would be if the second helix were ratio- 
nal. However, after many turns the angles at which the 
discrete charges on the two helices are placed will begin 
to drift with respect to the corresponding positions for 
a rational pair of helices. This is because the -0 value 
for the irrational case is very close but not equal to the 
rational tp value. It follows that the irrational interaction 
energy is the 9 average of the rational interaction energy. 



See Fig. 3(a) 



We can look at this averaging from another perspec- 
tive, keeping 9i fixed and allowing the other variables to 
adjust. It is easy to see that through this averaging 
will rotate through 27r radians while A4> and C will re- 
main fixed throughout. Sec Fig. |3(b)] It follows that the 
energy for an irrational pair is the average of a nearby 
rational pair with the same parameter values for Aip and 

C- 

If we take the Aip average of the irrational energy we 
obtain the interaction energy between one continuous he- 
lix and a second irrational helix. This is equivalent to the 
interaction energy between two continuous helices, how- 
ever, because the set {ki/j : /c G Z} is dense modulo 2tt 
for irrational ip. We may, therefore, smear out the charge 
over the first helix's backbone as well without changing 
the energy. It follows that the continuous pair interaction 
energy is the Aip average of the irrational pair energy. 

Finally, taking the ^ average of the continuous case, 
we obtain the interaction energy between a continuous 
helix and a cylinder with charge uniformly distributed 
over its surface. The potential due to the cylinder is 
independent of the position along the direction of its axis. 
The energy is thus unaffected when the continuous helix 
is averaged out in this direction, as well, and so this is 
equivalent to the interaction between two cylinders of 
charge. The cylindrical model interaction energy may, 
therefore, be obtained from a C average of the continuous 
system's energy. 

One significant consequence of the above follows from 
the fact that if one system's energy is the average of an- 
other's, the latter's energy can take on values both larger 
and smaller than that of the former's through an appro- 
priate choice of the parameter being averaged over. Phys- 
ically, both the irrational and rational situations may be 
realized. As the irrational case is the (j> average of the ra- 
tional energy, it follows that an interacting pair of helices 
can often decrease its interaction energy through an ad- 
justment of the angle ip to a nearby rational value. In the 
following section we demonstrate that the amplitude of 
the (j) dependence decreases exponentially with the value 
of . This statement specifies the manner in which the 
rational energy benefit depends on the rationality of ip. 
With this information we can now determine whether or 
not all irrational helices will be induced to twist when 
interacting with a second identical helix. 

We start by positing that the energy gain associated 
with a given rational helix pair scales as exp[— 771^]. 
Further, we suppose that this energy gain is sufficient 
to twist all irrational helix pairs in the twisting range 
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(a) 



(b) 



FIG. 3: 



3(a) Because of angular drift, the energy for an irrational pair will be given by the 6 average of a nearby rational pair. 

average, with A?/; and C, kept fixed. In 



|3(b) Taking the first helix's perspective of the 6 average, we see this is equivalent to a 
the figure 



(-0 — Sip^ tp + Sip) to the rational value tj:. It follows that 
dip is also proportional to exp[— 771.0]. If we sum up the 
widths of the twisting regions around each rational tp 
value less than 27r, we get 

= ^ exp[-7n] 

m/ne{Q<l} 

00 n 

< ^ ^ exp[-7n] 



n— 1 m— 1 



(e^ - If 



(6) 



In the second line we have replaced the sum over all ra- 
tional values by the sum over all relevant integer denomi- 
nators and numerators. Note that this significantly over- 
counts the number of rationals since we are including 
terms which are not in reduced form. Additionally, wc 
overestimate the portion of the interval associated with 
rational helices, in that some of the basins of attraction 
of rational helices are "shadowed" by others. In spite 
of this over-counting, the sum over the twisting widths is 
bounded. Indeed, for large separations, we expect 7 to be 
large. In this case VF(7) < 27r, implying that not all the 
irrational ■0 values in (0,27r) will be twisted. However, 
as the helices approach one another, more and more irra- 
tional helices should experience a (p induced twist since 
the energy benefit of a rational ip) value increases with 
proximity. This expectation is confirmed when one plots 
the optimal angle ip for an interacting pair of helices ver- 
sus the isolated -ip value at different spacings d. This is 
done in Figs. 4(a) and 4(b) At large d one observes an 



incomplete devil's staircase. As the separation distance 
is decreased, a filled in staircase is observed and nearly 
all helices are twisted to a nearby rational value. 

A second significant consequence of the above averag- 
ing arguments is that they often allow for a convenient 
decomposition of the interaction energy, as mentioned 
above. Beginning with the cylindrical model interaction 
energy, we may add on corrections for the continuous he- 
lix, discrete irrational, and finally discrete rational terms. 



Each new correction term adds dependence to the energy 
on a new parameter. In addition, at the level at which 
energy dependence on a given parameter is first intro- 
duced, the energy has a well defined periodicity in that 
parameter. For example, although the irrational energy 
is not periodic in the continuous energy is and has the 
finite period of 27r. The periodicity at each level allows 
each of the correction terms to be expanded in a Fourier 
series. Further, except in certain extreme limits, these 
Fourier expansions are quickly damped. Truncation of 
these series thus allows for simple, approximate expres- 
sions for the energy to be obtained which are consistent 
with the sinusoidal forms typically observed numerically. 
While first order truncations are often sufficient, higher 
order harmonics may be required to accurately model 
the energy for parameter locations which allow for close 
charge interactions. This is often the case for C, values 
near tt and for small d values, for example. In this case 
the helical backbones are close to one another and small 
adjustments in l^ip or (p may allow for large increases in 
the energy. In practice, physical systems should often be 
found far from such parameter locations, however, and 
first order approximations for the energy are therefore 
acceptable. 

Numerical observations have indicated that the contin- 
uous energy often dominates both the A-0 and (p depen- 
dence by at least one order of magnitude. In this case, 
we say that the system is in the continuous limit and 
think of the rational and irrational terms as perturba- 
tive correction terms for the energy. For small values of 
the parameter d, however, the correction terms can also 
have amplitudes on the order of magnitude of the ther- 
mal ksT energy scale per azimuthal persistence length 
[l0| . Notable extreme limits where the continuous limit 
is not valid include the large a limit, in which the charges 
on each helix are separated by large axial distances, and 
the small a, rational limit, where each helix looks much 
like a grouping of lines of continuous charge. Most 
physical systems appear to be somewhere between these 
two limits, however, and the energy may be considered 
to be in the continuous limit with first or second order 
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sinusoidal discrete correction terms sufficient. 



IV. PHASES AND AMPLITUDES 

To continue the characterization of the pair interaction 
energy, we now consider the phases and ampUtudes of 
the Fourier expansions of the correction terms described 
above. In addition, we briefly discuss the irrational en- 
ergy's "0 dependence. This "0 dependence is of interest 
since it may provide a mechanism for discrete interac- 
tion induced twisting for some systems. 



Continuous C, dependence 

The phase of the continuous C dependence may be de- 
termined by explicitly writing down the energy in integral 
form and differentiating to find extrema. Doing this we 
find 



Here, 



and 



2?- 



A(2^Ph^)^d^,d^,. (7) 
dx X dQ 



dx dx dx"^ 
dC " d^~d^ 

cx a^ie^-e^ + O. (8) 



a^{92 - 6*1 + C)^ - 2rd(cos6li - cos 6*2) 

-2r^ cos(6li - 02).{9) 



The integrand is odd about Oi = ( and 62 = when ( = 
kir, with k an integer. The continuous energy correction 
term may therefore be expanded as 



(10) 



k=l 



The magnitudes of the coefficients A^x have been con- 
sidered in previous studies of the continuous energy 

and have been proven to decay exponentially with 

both d and k [1^. The decay rate was found to be 
exp[— 27r-\/fc2 -I- a-^d/P], where P is the pitch of the he- 
lices. We note that this result could also be obtained 
through a superposition of interactions between lines of 
periodically spaced point charges. The rapid convergence 
of the above sum leads us to expect the first coefficient 
Ai^i^ to take on a negative value quite generally. 



Irrational Ail) dependence 

Explicit integral expressions for the Fourier coefficients 
of the irrational energy may be obtained as follows. Re- 
calling that the set {kip : k G Z} is dense modulo 2tt for 



irrational ip, we set rii = n2 + ma in Eq. The sum 
on 712 rnay then be replaced by an integral over to 27r. 
Next we apply the Poisson sum rule to the sum on rria 
and obtain for the energy per charge on the second helix 



Er 



q 



2ttiP J_ 



27r ^ — a_^x 



-deTdO: 



(11) 

where x^ is again given by Eq. Note that the fc = 
term above gives the continuous energy. Upon averaging 
over Alp the other terms vanish, an observation consis- 
tent with the above geometric averaging result. 

Parity arguments again allow us to determine the ex- 
trema of the irrational Aip dependence, but only when 
C is an integer multiple of tt. At integer C/tt, it may be 
easily shown from Eq. (|lip that there are irrational Ai/j 
extrema at Aip = C, C ± ■0/2, C ± 20/2, etc. Although 
symmetry arguments alone are insufficient to determine 
the phases away from integer C/""; numerical observa- 
tions indicate that this phase is often quite linear in C. 
Sec Fig. [5l The slope of this phase dependence in ( de- 
pends upon the parameters a, r, and d, however. This 
information allows us to expand the irrational energy as 



27rfc 

Ea^, = ^ Ak,A^ cos(^p[A0 ■ 

k=l ^ 



(1 



ZTT 



■o(C')]), 



(12) 

terms in the phase 



where s is some integer and the 0(C^ 
must vanish whenever i^/tt is integer. 

We may demonstrate that the Ak.A^i in Eq. (fT2|) decay 
exponentially with k by returning to the integral repre- 
sentation in Eq. (|lip . Changing variables to (u, w) = 
{62 + 61,62 — Oi), it may be shown that for d > 2r, one 
may always add to w a positive, finite imaginary part 
iK without crossing any singularities. For large d, n is 
bounded by d/a and 

. . 2'Kkd^ , , 

Ak,A-4i oc exp — , (13) 

aip 

also consistent with the discrete lines of charge limit. 



Rational dependence 

Once again setting ni = n2+ma in Eq. the rational 
energy per charge on the second helix may be expressed 
as 



where 



^ nia — — oc k—1 



d^ + 2r^ + a^{maip - Aip + Cf 



(14) 



-\-^rd sin( 



ruaip + Alp 



27rfc . ruaip - Alp 
)sin( ) 



-2r^ cos{maip — A0). 



(15) 





(a) 



(b) 



FIG. 4: Shown are two plots of the preferred angle between charges ^ for an interacting pair of helices versus the value of 
^ assumed for an isolated helix. The helices were given an arbitrary torsional spring constant which resists twisting, which 
is assumed to occur in such a way that the linear charge densities of the helices is fixed. In 4(a) the separation distance is 
large enough so that only an incomplete devil's staircase is observed. In 4(b)[ the separation has been decreased and a filled 
in staircase is observed. Note that the one charge per turn system was excluded. If this is retained, it dominates the small ip 
region of the plot and all irrational helices are observed to twist. 
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FIG. 5: A contour plot of the irrational energy versus A?/> and 



It is easy to sec that the above is even in ((> about the 
point (At/i — 7r)/2, independent of C,. To show this, one 
need only note that vria'^ll = TOaTi^Tr/n^. This phase 
shift will not affect the parity of inside sum over fc since 
it always shifts the sum by an integer multiple of 7r/n^. 
It follows that the rational energy correction term may 
be expressed as 



= ^ ^fe,</> cos[n^fc((^ - 



)]■ (16) 



To examine the convergence of the sum in Eq. p6p . 



we again focus on the inside sum of Eq. ([T^ . Letting 
(5 = 0— "lai/M-Ai/j ^ rewrite this inside sum as 



/("^» = + ^cos( 



27rfc 



S)). 



(17) 



Here a and b are constants over the k sum, with b < a 
since d > 2r, and g is the individual charge potential. To 
obtain an expression for the first coefficient in Eq. (jl6p . 
we Taylor expand the function g about the point a 



j 



(18) 



and note that only those j > In^ contribute to A;,^. 
This is because to get an argument containing ln^6, you 
need a product containing at least In^ factors of cos(<5). 
Plugging in the Coulomb potential and summing on j > 
In^ gives 
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y cos(— - 6)'"^+' oc i-Y"^. 



(19) 

It may be shown that the fluctuating portion in the co- 
sine sum in Eq. ()19p is of order unity. Since 5 < a, it 
follows that the series Eq. (fTB|) is exponentially conver- 
gent and that the decay rate is proportional to n^. A 
careful look at the values of b and a above shows that 
the (j) dependence decays algebraically with d. Thus, at 
large distances, the rational correction terms will dom- 
inate the continuous helix correction terms. Typically, 
this will occur at sufficiently large distances that each of 
these corrections are effectively negligible. For small a 
this limit may become significant, however. 
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(a) 



(b) 



FIG. 6: 6(a) A plot of the rational energy versus (/> and A-i/) for -0 = 27r/7 and C — 0. Note that the first component vanishes 
when Ai[) is an integer multiple of ■(/), consistent with Table |l] |6(b)| A plot of the rational energy versus (p and for ip = 2tv/5 
and Atf) = 0. While the phase of the (p dependence is independent of (, the amplitude is not and is observed to change signs 
as ^ moves through 2tt. 



It turns out that the phase dependence in Eq. (fT6|) may 
sometimes have an interesting impact on the amplitude 
of the coefficients. Consider what happens when (jj and 
Atp are each adjusted by one period. Plugging into the 
k = 1 term in Eq. (fT6|) gives 

A^' + - n 
^i,4AV + ^, C) cos[n^((/) + ^ )] 

= (-1)™* Ai,4A^ + ^, cosK(</) - ^^^)]. 

(20) 

From periodicity, the energy must be unaffected by this 
shift in (f) and Aip. It follows that 



A,,4,{AiP + ij) - (-i)'"*^i,4A7/.). 



(21) 



Therefore, for m.^ odd, Ai^^ must change sign as Aip 
is adjusted through -0 radians. Since Ai_^ is a smooth 
function in Aijj, it follows that for odd there is a Atp 
value at which the amplitude of first component van- 
ishes. This is significant because the second component 
is in general exponentially smaller than the first compo- 
nent. The result is that there is little dependence at 
this value of A-0. See Figs. 6(a)| and |6(b)[ 

In general, the particular A-ip location where the first 
component vanishes depends on the values taken by the 
other parameters. For certain high symmetry ( values, 
however, the exact A'0 values may be determined and are 
independent of the remaining parameters. To determine 
these locations, wc must note that the energy is invariant 
under the operation 



TABLE I: Solutions to Eq. (|23|l which specify exact param- 
eter locations where first component's amplitude is zero. 
Here, the ki are arbitrary integers. 





n.4, 


C A^A 


odd 


odd 


± kinm^ ± k2ip 


odd 


even 


± fcivrm^ ^ ± k2ip 


even 


odd 


± fciTrm^ ± k2tp or ^ ± k^ip 



(0, A0, C) ^ (0 + TT - AV-, -AV-, -C). 



(22) 



This operation is equivalent to relabeling the helices one 
and two. If the first component vanishes, this suggests 
that wc look for parameter locations where the period 
of the dependence is halved. That is, we search for 
parameter locations where the energy is unaffected when 
is adjusted by tt/ti^. Plugging into Eq. ([^ we look 
for solutions to 

E{<Po,Ai;,0 = £;(0o + 7r/n^,A0-,C) (23) 
= £'(0o + 7r/n^-f^7r- A'0, -A'0, -0- 

The solutions to Eq. ((23)) are shown in TableU Note that 
even solutions exist as well. While the amplitude 
does go to zero at these even locations, it does not 
change signs. Further, while the amplitude of the first 
component is required to vanish at other ^ values for odd 
, the even solutions noted in Table [T] are the only 
solutions observed numerically. This is consistent with 
our expectations: while the zero-frequency component of 
the Fourier expansion of must vanish for odd m^, 
symmetry conditions do not require this to be the case for 
even m^. The expansion of for odd is discussed 
later in our consideration of A-DNA. 
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Irrational ^ dependence 

Due to the properties of the rational ip dependence, 
it is clear that the pair interaction energy is a nowhere 
continuous function in ■0. In addition to rational energy 
terms there is also a direct irrational ip dependence. This 
irrational ip dependence provides another possible mecha- 
nism which may be responsible for the observed twistings 
of helical macromolecules in aggregate. 

The functional form of the ip dependence will depend 
on how the helices twist or untwist. In the case of F-actin, 
the heights between charges remain roughly fixed during 
twisting [2^ . This implies that the linear charge density 
of the molecules is unaffected by the twisting. Numerical 
plots demonstrate that for systems which twist in this 
way, the locations of the irrational energy minima in ip 
depend strongly on the parameters d/r and a/r. There- 
fore, to determine whether this mechanism is a possible 
cause for twisting, one must numerically examine the ap- 
propriate phase space region for the system of interest. 

Resulting energy expression 

The general energy expression, then, is given by adding 
to the energy of two interacting cylinders of charge the 
expressions in Eqs. (fTO|) . (fT2|) . and (fT6)l . It is important 
to remember that the A^,^ depend on all of the param- 
eters except (j), the Ak_Atii depend on all the parameters 
except for (p and Aip^ and that the Akx depend on all 
parameters except for (p, Aip, and The rule of thumb 
is that a particular amplitude will be large when close in- 
teractions may be introduced through the particular pa- 
rameter's adjustment. Although we focused on locations 
where the first component of the (p dependence vanishes, 
the Alp dependence amplitudes can also take on both 
positive and negative values depending on the values of 
the parameters a and d. In general, we must resort to 
numerics to determine which terms dominate the energy 
landscape for a given system and what the signs are for 
the amplitudes of the various terms. 

V. MODE ANALYSIS 

We turn now to a consideration of helices which are not 
rigid but instead have some internal degrees of freedom. 
We first consider the modes of a system of mobile charged 
particles constrained to move on the surface of an isolated 
cylinder. Each particle has the same charge and the axis 
of the cylinder contains a compensating line of charge of 
the opposite sign which ensures a net charge neutrality 
for the system. We assume an initial helix distribution 
for the charges, which by symmetry, is clearly stable with 
respect to the motions of any single charge in the system. 

For simplicity, we consider modes in which the charges 
are only allowed to rotate in the azimuthal direction. 
Thus, the axial positions of the charges arc fixed. Recall 



that this is roughly how the charges fluctuate in F-actin 
systems. For a given mode, the position of the charges 
can be written as 

f„ = Rn+u„, (24) 

where i?„ describes the equilibrium position of the charge 
n. This is given by 

Rn ~ r(x cos nip + y sinnip) + zhn, (25) 

where we have introduced the notation h = aijj for the 
axial rise per charge. The vector Un is the displacement 
from equilibrium and is given to second order in the mode 
amplitude Aq by 

M„ = {Aqe'""'' + A*er"^''''){-xsinniP + ycosniP) 

— {xcosmp + ysv[inip).{2&) 

We assume a Coulomb interaction between charges and 
make use of the identity 

/•oo 

= 7r-i/2 / u;-i/2e-«"-'rfu;. (27) 

We desire an expression for the energy valid to second 
order in Aq. To that end we Taylor expand Eq. ([TT)) 
with respect to r up to second order and then sum up 
the contributions from each pair of charges n and m. This 
is a straightforward but lengthy procedure and we quote 
only the result. Throwing out negligible terms which do 
not scale with the number of charges in the system, we 
obtain the change in potential energy 

~ ^ —4(1 — cosqh{n — to)) cos(n — m)Tp 

n.ni 




(28) 

Note that the resulting expression depends on the indices 
n and to only through their difference. This means that 
a sum over n with n — m fixed yields a multiplicative 
factor going as the number of charges in the helix. Thus, 
the potential energy per charge is given by summing over 
k = n — m in Eq. (|28p . This sum was re-expressed in a 
form suitable for quick numerical evaluations using the 
Ewald summation approach, just as in the helix-helix 
interaction calculation shown in the appendix. The pro- 
cedure is again straightforward. However, the resulting 
expressions are lengthy and provide little intuition and 
so will not be presented here. Two plots of the resulting 
dispersion relations are given in Fig. [71 In each case the 
initial angle between charges was chosen to be V' = i". In 
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the first case h = 3. Or and in the second case h = 1.8r. 
What we see in the first plot looks like a typical phonon 
dispersion relation. As the wavevector increases, the fre- 
quency increases from zero. In the second plot, however, 
we see that there is no real frequency solution for small 
wavevectors. This indicates that the = t: system is 
unstable with respect to global rearrangements of the 
charges for this value of h. 

The instability of the ip = tt system at small h/r val- 
ues can be understood through a consideration of the en- 
ergy of a single helix of charge. Because the self-energy 
becomes variable for flexible systems, this energy must 
also be taken into account when considering inter-helical 
interactions. In order to calculate this energy for an iso- 
lated, perfect helix, we once again carry out an Ewald 
summation. We sum up the contributions from each 
charge to the potential at a given point on the cylin- 
der. We then take the limit as this point approaches the 
position of one of the charges and subtract off the inter- 
action with that charge to get the potential due to each of 
the other charges at this location. The divergence of the 
potential due to the infinite number of other charges on 
the cylinder is canceled out when added to the potential 
from the neutralizing charge distribution up the center 
of the helix. The resulting expression for the energy per 
charge is 



what the charges would look like if we were to unroll the 
cylinder. What we would see is what appears in the rect- 



ht 



ht 



l,m' 



4E r^'^ h{2rH)er'^''-'^'T^t-'dt. (29) 







In the above, /3 is the constant determining the cut-off 
between the high and low integration domains, the /; are 
Bessel functions, and the primed sum is over all I and m, 
excluding the term where they are both zero. 

Resulting energy versus plots for a single helix are 
shown in Fig. [51 At large values of h/r a single energy 
minimum appears at two charges per turn. However, as 
the value oi h/r decreases below approximately 2.1, a bi- 
furcation occurs and two new minima replace the original 
minimum at ip = n. A slight twist is introduced, either 
right or left-handed, in order to increase the separation 
distance between axially adjacent charges on each side of 
the cylinder. This explains the observed instability of the 
Ip = TT system at small h/r. As h/r is further decreased, 
more and more nearly degenerate minima appear. This 
is explained below. 

At small h/r, the charges are tightly packed onto the 
surface of the cylinder. They will thus attempt to ar- 
range themselves into a structure resembling an energy 
minimizing triangular lattice. To determine which ip val- 
ues allow for nearly triangular lattice packings, consider 



angle in Fig. 9(a) The rectangle is a portion of the rolled 
out cylinder, which continues to the left and right. The 
charges on the cylinder, viewed in this way, are a portion 
of a Bravais lattice of charges in which the spacing be- 
tween vertical lines along which the charges lie is h and 
the distance between charges on one of the lines is 27rr, as 
shown in the figure. We may interpret the lines on which 
the charges lie as Miller lines of the lattice. Note that the 
volume of the primitive cell of this lattice, vws — ^irrh, 
is completely determined by r and h. Consider now a 
triangular lattice with this same primitive cell volume. 
Suppose this triangular lattice has a Miller line on which 
the charges are spaced by 27rr. This would determine a 
helical configuration for our cylinder of radius r and ax- 
ial rise h which would unroll into this triangular lattice 
structure. The angle ip between charges would be deter- 
mined by the spacing between charges on adjacent Miller 
lines. In general, there will probably be no Miller line 
on the triangular lattice with charge spacing 27rr. How- 
ever, there may be Miller lines with spacings which are 
close to 27rr. If this is the case, a slight distortion of the 
triangular lattice would result in a possible structure for 
the helix which is energetically favorable. To determine 
candidate values, one can take a point on the triangular 
lattice and draw a circle of radius 27rr about this point, 
as shown in Fig. 



9(b) 



At small h/r, the radius will be 
large compared to the charge spacing and many points 
will appear near the circle's outer perimeter. Such points 
provide candidate Miller lines which upon a slight distor- 
tion will determine a helical configuration of charges that 
results in a nearly triangular lattice. These indeed corre- 
spond to the energy minimizing ip values. An example is 
shown in Fig. [TUl We can thus understand both the small 
and large h/r limits of the single helix energy landscape. 

We now consider the interaction between a pair of flexi- 
ble helices. Once again we assume that the axial positions 
of the charges are fixed but allow them to independently 
rotate in the azimuthal direction. All charges interact via 
a Coulomb force law and all charge values are taken to 
be 1. In addition, an clastic energy is introduced which 
resists the adjustment of the azimuthal separations be- 
tween nearest neighbor charges on each helix. The en- 
ergy cost of adjusting a nearest neighbor separation is 
kA9^/2. Here, k is an elastic constant and A6 is the 
difference between the assumed angular separation for 
a nearest neighbor pair and ip, the equilibrium angular 
separation assumed when the helices are isolated. 

In order to determine the mode structure of a given 
pair of helices, one must first find the equilibrium ori- 
entation. This equilibrium orientation is characterized 
by a given set of global parameters, but also requires 
a determination of the optimal internal twistings within 
each helix. We employ the matrix version of the Newton- 
Raphson method to determine these equilibria. We start 
by assuming a given orientation for the two helices. This 
is specified by the parameters (f>, Aip, and ip, where ip 
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FIG. 7: Dispersion relation plots for the helix modes of two = n systems. In the first h = 3. Or and in the second h — l.Sr. 
In this second case, the frequencies for certain long wavelength modes were found to be complex indicating an instability. 
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FIG. 8: Plots of the single helix energy as a function of ij] for the various values of h/r indicated in the figures. 
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FIG. 9; 9(a) The spiral arrangement of charges on the rolled out cylinder. 9(b) Candidates for the Miller line 
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2.83595 
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FIG. 10: Shown are the energy minimum a.t tp — 2.83595 which appears at h/r — 0.05 and the rolled out lattice structure it 
corresponds to. 



is again taken to be 2Trm^ /n^ . Periodicity is enforced 
every charges for each helix. Thus, once the distance 
between the helices and the axial positions of the charges 
are set, the geometry of the system is determined by the 
vector 6, which contains 2n^ components specifying the 
angular positions of the charges on the two helices. The 
force vector and the Hessian matrix for the system is then 
calculated. Their components are defined as 



m 



dE 

89,89, ' 



(30) 
(31) 



where E is the total electrostatic and elastic energy for 
the system. If the force vector does not vanish identically, 
the angular positions of the charges are adjusted by a 
small amount 59 in order to reduce this force. Assuming 
the system is near an equilibrium point, one can take a 
first order Taylor series for the force and set it equal to 
zero to obtain 



orientations, the amplitudes of the motions of the two 
helices are equal for each of these modes. In this case, 
the first mode corresponds to adjusting only the (f) global 
degree of freedom. From our symmetry arguments above, 
we know that the interaction energy of an irrational pair 
of rigid helices does not depend (/>. We thus expect there 
to be a gap between the irrational and rational frequen- 
cies of oscillation in this mode. In order to test whether 
or not this could be observed in flexible systems we com- 
puted the oscillation frequencies of this mode for various 
integer number of charges per turn systems. The elastic 
constant was taken here to be rather large such that little 
internal motion would occur. A log plot of the resulting 
oscillation frequency versus n,p is shown in Fig. 1111 Note 
that the log plot is approximately linear, consistent with 
the exponential decay in of the 4> dependence for rigid 
helices. We conclude that for systems with low flexibility, 
a gap should be observed between irrational and rational 
helices at long wavelength. 



H,^,59j - /. « 0. (32) 

An approximation to the ideal 56 is then obtained after 
a matrix inversion. This is added to 9 to obtain a new 
9 and the process is repeated until it converges to an 
equilibrium location. 

In order to determine the mode structure of the sys- 
tem one need only calculate the eigenvectors and eigen- 
values of the Hessian matrix of the energy. Note that 
this is conveniently already determined after iteration of 
the Newton-Raphson method. Of particular interest are 
the two lowest energy modes of the system. These modes 
involve the rotations of the two helices about their respec- 
tive axes with little internal degrees of freedom excited. 
Indeed, in the rigid limit, these are the only two modes 
allowed since the relative axial positions of the charges 
are held fixed. In one of these modes the helices rotate 
in the same azimuthal direction while in the other they 
rotate in opposite directions. For symmetric equilibrium 



As the flexibility of these systems is increased local 
twisting will start to become more significant. The two 
helices may then be modelled as a pair of perfect helices 
plus a series of physical dipoles. If the dipole strengths 
are weak, or if the system is in an orientation where the 
rigid helix azimuthal energy dependence is large, the per- 
fect helix interaction will dominate. However, as the flex- 
ibility is further increased, dipolc-monopole interactions 
can begin to dominate the azimuthal energy dependence. 
An example is provided by the seven charge per turn sys- 
tem depicted in Fig. [T^l At and below a threshold value 
of K, the stable equilibrium orientation switches to the 
second orientation shown. This corresponds to a change 
in the equilibrium (j) value of 2it /I . This example demon- 
strates that local twisting can alter the interaction in 
a qualitatively considerable way even for systems which 
contain only a modest degree of flexibility. 
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(j), Alp and C- The effective energy expressions will thus 
J allow us to model many physical aggregate systems with- 

• ^ out having to know the detailed form of the interaction. 

0.1 There is one caveat, however; the interaction energy in a 

system of flexible helices is no longer pair additive, due to 
• the fact that helices may interact indirectly through the 

_ local twisting of third party helices. In the applications 

• below, wo assume that the pair interaction dominates 

r* • such effects. This assumption may not hold for highly 

flexible systems. 



It 



FIG. 11: A plot of the oscillation frequency of the sym- 
metric mode versus n^. In this mode the helices rotate 
about their respective axes in the same direction. Here, 
the system parameters were defined by d = 2.5r, = 1, 
a — l/^f) — 2n /n^, and C = 0. 



VI. PHYSICAL APPLICATIONS: F-ACTIN AND 
A-DNA IN AGGREGATE 

Both F-actin and A-DNA have been observed to con- 
dense under the influence of multivalent counterions. 
These two systems occupy different limits in the energy 
landscape and provide relevant systems to which we may 
apply our effective energy expressions above. In deriv- 
ing these expressions, we assumed that the interacting 
charges formed perfect helical distributions. However, 
these same expressions also often apply for disordered, 
flexible systems. To prove this, one must note that the 
pair interaction is invariant under the operation 

(</),A^,C)-^(-0,-A^,-C)- (33) 

This operation is equivalent to rotating the helices by 
180 degrees so that they are flipped upside down. This 
symmetry, together with that in Eq. ((22|) . allows for an 
immediate, and more general, derivation of the phase re- 
lationships observed in Eqs. HH), and ([ill). This 
proof applies whenever the variables (j), Atp and C a-rc 
well-defined and provide a sufficient characterization of 
a pair's mutual orientation. In particular, it applies if 
there is uncorrelated, non-additive disorder in the charge 
locations and it also applies if the system is flexible. Note 
that flexibility introduces many more degrees of freedom 
to a system. However, there should remain a single well- 
defined ground state orientation for any choice of the 
three parameters (j), Aip and C, and it is the energy of 
this ground state which will be represented by the effec- 
tive energy expressions considered here. Note that this 
is consistent with the change in equilibrium orientation 
observed in the flexible system depicted in Fig. [T^l ap- 
parently increasing the flexibility in a system can allow 
for the amplitudes of the various terms to change signs, 
but the phases must remain fixed. 

The argument outlined above demonstrates that the 
effective energy expressions derived in section IV apply 
whenever the pair interaction maintains the periodicities 



A. F-actin 

The average ip value in Mg^"*" condensed paracrystals 
of F-actin is 27r6/13 w 2.9 radians. This value differs 
from the average ip value observed for isolated actin fila- 
ments by about 0.5 degrees per monomer j^^]. Numerics 
indicate that for a discrete ip = 27r6/13 helix pair, with 
radii and a values scaled appropriately for actin, the con- 
tinuous energy dominates globally. That is, the energy 
is minimized at, or very near C = 0. Further, at = 
it is found that the irrational energy terms dominate the 
rational energy terms and the extremum at Atp = is 
the location of the energy minimum. This is consistent 
with the observation that filaments within a given layer 
of these paracrystals are typically found to be in register 
[2l|. Taking ( = Aip = to be the location adopted by 
the physical system, it follows from Table U that the am- 
plitude of the first <j) component correction term is zero. 
This will have the effect of greatly reducing the energy's 
dependence on the parameter </>. 

Numerics also suggest that the irrational ■0 dependence 
cannot alone be responsible for the observed twisting of 
F-actin in aggregate. For small d, the energy minimum 
occurs at = TT, which would always tend to increase "0. 
This is inconsistent with the observation that different 
types of actin were each observed to twist ioip = 27r6/13, 
regardless of whether the isolated -0 values were less than 
or greater than this value At larger d values, near 
d K, 5r, a moving energy minimum appears near this 
value. At these distances the energy gain associated with 
the irrational energy is insigniflcant compared to the en- 
ergy cost due to twisting, which was evaluated using the 
torsional rigidity constant reported in Ref. [2^. It thus 
appears unlikely that charge discreteness effects can al- 
low for sufficient energy gains to explain the 0.5 degree 
twist per monomer observed in Ref. (20| . 

More recently, tightly packed Ba^"*" condensed F-actin 
bundles have been observed which have been torsionally 
twisted to an average ip value of 27rl7/36 [1^. Note that 
the results of Table |T] indicate that the first (p compo- 
nent would not vanish for a C = A?/; = hcxagonally 
packed lattice of molecules of this configuration. The 
value n.^ = 36 suggests that the rational correction terms 
should be highly damped for this system, however. Nev- 
ertheless, without knowing the precise form of the in- 
teraction wc cannot rule out that this twisting was in- 
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FIG. 12: Shown are end-on views of two equilibrium orientations for a seven charge per turn pair. The larger black dots 
represent the locations where the charges would sit for a rigid system at equilibrium and the smaller blue dots represent the 
corresponding locations for the flexible system. For a rigid system, the first orientation would be stable. However, for elastic 
constants below k = 9, the value for the system depicted here, the first orientation becomes unstable due to local twisting and 
the second orientation shown becomes the stable equilibrium. Here d = 2.2r, -0 = 2n/7, a — l/'0, and C = 0. Note that the two 
charges which appear to be very close to each other in the second figure are actually separated by a significant axial distance. 
(Color online) 



duccd in order to obtain a rational energy benefit. Actin 
monomers are highly heterogeneous, and as the authors 
of Rcf. [2^ point out, their highly-charged subdomain-ls 
may often dominate the interactions. 

It is important to point out that F-actin is particularly 
flexible in the azimuthal direction. It is thus possible 
that the inclusion of non-pairwise additive local twist- 
ing effects might help to explain the observed averaged 
deformation angles. 

B. A-DNA 

Previous numerical work has shown that for both 
A-DNA and B-DNA the correction terms due to dis- 
creteness cannot be mutually optimized for each nearest 
neighbor pair in a hexagonally packed lattice [l^ . Due to 
frustrations in the C dependence, B-DNA is typically ob- 
served to pack in an orthorhombic lattice, however [2^ . 
The resulting reduction in the number of nearest neigh- 
bors for each molecule makes it more likely that frustra- 
tions in the rational dependence can also be overcome. 
As discussed below, there is no C frustration in aggre- 
gates of A-DNA. However, A-DNA is also not observed 
to form hexagonally packed aggregates [l^. In order to 
strengthen the plausibility of the suggestion in Ref. 
that other lattice structures are assumed in order to re- 
duce azimuthal frustration, we shall now revisit the prob- 
lem of packing A-DNA on a hexagonal lattice. We find 
that to lowest order the frustration results in a complete 
averaging out of the dominating azimuthal, discrete en- 
ergy terms. 

Because A-DNA is a double helix, we must consider 
four single-helix pair interactions for each pair of neigh- 
boring molecules. Each of the single helices of A-DNA 
is observed to have 11 charges per turn in aggregate. 
Experiments and numerics indicate that the energy is 
optimized when there is no axial shift between the two 
molecules [l3|- This implies that two of the single-helix 
pairs will correspond to C = interactions. The remain- 



ing two interactions, which correspond to the upper he- 
lix of one molecule interacting with the lower helix of 
the other, will have non-zero C values. At all ^ values 
the rational energy terms are observed to dominate the 
irrational energy terms and we may approximate the az- 
imuthal energy dependence by the first order rational 
correction to the energy. Because the energy is periodic 
in At/i with period ip, we may expand the amplitude of 
the first (j) component, ^1,0, as 

^ . /ll(2fc+l)A7A \ 
^1,0 = "fc ( ^ + 7fc J ■ (34) 

k=Q ^ ^ 

Numerics indicate that the first term alone often accu- 
rately models the interaction and we have the follow- 
ing approximate expression for the pair interaction's az- 
imuthal dependence. 

,11 Alb , , n-Ath, , , 
S^ = asin(^-^+7)cosll((/)+ ^). (35) 

Here, both a and 7 are functions of C. 

Consider first the interactions between the upper helix 
strands. For these interactions C = and the results of 
Table [J indicate that 7 = 0. Plug ging in Eq. (|35p for 
each of the six nearest-neighbor interactions for a given 
molecule demonstrates that the energy is independent of 
that molecule's azimuthal orientation. This implies that 
these terms average out in the bulk. For the same reason 
the interactions between the lower helix strands must also 
average out. 

Now consider the two non-zero C interactions for a 
given pair of molecules. If one of the interactions cor- 
responds to = (ijD the other must correspond to 
C = Cdu — —C,uD- Here the subscript UD indicates the 
interaction between the upper strand of the first molecule 
and lower strand of the second, while the subscript DU 
indicates the opposite interaction. To relate the a and 
7 values for these two interactions we apply the trans- 
formation appearing in Eq. (j22[) which leaves the energy 
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invariant. 

a{CuD) sin[— h liCuD)] cos 11(0 H ) 



= a(CDC/)sin 



.llAi/) 



TT — Atjj 

- 7(C(7n)]cosll((/iH 



It follows that 



aiCuo) 
liCuo) 



~i{Cdu) 



2 

(36) 



(37) 
(38) 



One final geometric effect must be taken into account. 
In general, the azimuthal positions of the charges on the 
upper strands may be shifted with respect to the po- 
sitions on the lower strands. Let the mean shift be 5 
so that for two corresponding charges on the upper and 
lower strands of one molecule we have, 9i_u = Oi^o + ^■ 
With this definition the sum of the two non-zero C, inter- 
action terms is 



Ed: 



r . , ll{A^J-S) , . , ll(AV> + ^) , 
"i sm[ h7j +sm[ 7j 



cos 11( 



TT — Alp — S 



(39) 



= 2acos(-^ 7) sin — - — cos 11(0 H )(40) 

where = 0—^/2. This has the same form as the C = 
interactions and, therefore, also averages out in the bulk. 

In the Kornyshev-Leikin theory of helix-helix attrac- 
tion, a large number of counterion charges are assumed 
to bind to the grooves of the helices [l3|- This groove 
binding allows for the correlations necessary for attrac- 
tion. To take any groove bound charges into account, we 
assume, as in Ref. [l^, that the groove bound charges 
arc not azimuthally ordered, and so taken together may 
be modelled as continuous helices of charge. The result- 
ing lowest order expression for the interaction energy be- 
tween one phosphate strand of one A-DNA molecule and 
a continuous, condensed counterion helix on a second A- 
DNA molecule takes the form ^ cos 11(0 -I- S). This en- 
ergy form averages out when summed over the six nearest 
neighbors of each molecule in a hexagonally packed lat- 
tice. Thus, taking all of the interactions into account, 
we have seen that to lowest order, the azimuthal energy 
terms completely average out in a hexagonally packed A- 
DNA system. This provides a strong statement regarding 
the degree of frustration for this geometry and is consis- 
tent with the idea that non-hexagonal packing structures 
are adopted in order to overcome azimuthal frustration. 



VII. DISCUSSION 

The work presented here is complimentary to that pre- 
sented in Refs. [l^l and [l^- The model considered in 



these previous articles attempted to incorporate the ef- 
fects of counterions explicitly. Here, we have focused 
directly on the symmetries of the interaction and have 
obtained some results which are, in a sense, model inde- 
pendent. In particular, the effective energy expressions 
derived can be applied to model both electrostatic and 
non-clcctrostatic aspects of the pair interaction. 

In our characterization of the helix-helix interaction we 
began by decomposing the energy into a sum of terms, 
each of which added dependence to the energy on a new 
parameter. For Coulomb interactions, the rational en- 
ergy correction term was shown to decay exponentially 
with n^. This result can also be shown to hold for any 
power law or Yukawa interaction. This characterizes how 
rational a helix pair has to be in order to obtain a signif- 
icant rational energy benefit. Symmetry arguments al- 
lowed us to demonstrate the existence of locations where 
the dominating rational correction term vanishes. The 
exact location in parameter space where this occurs was 
determined for certain high symmetry orientations of the 
two helices. Finally, the phases of the Fourier series 
expansions were determined for each of the correction 
terms. 

We next considered the mode and energy structure of 
an isolated flexible helix of charge. When the axial shift 
per charge is greater than h « 2.1, the single helix en- 
ergy is minimized when there are two charges per turn. 
Slightly below h « 2.1, axially adjacent charge interac- 
tions cause the energy to be minimized when the helix is 
slightly twisted either to the right or to the left. A series 
of similar energy bifurcations were observed to occur as 
h was further decreased. It follows that the energy of a 
single helix can be highly dependent upon twist angle. 
This may often play a role in determining the equilib- 
rium conformation of helical molecules. For interacting 
pairs of helices, the independence of the irrational en- 
ergy suggests that a gap should occur between rational 
and irrational systems at long wavelength. Scattering 
experiments could thus, in principle, provide informa- 
tion regarding the degree of rationality in a system of 
bundled helices. Although the phases of the interaction 
terms remain the same for flexible systems, it was shown 
numerically that the amplitudes can change sign. This 
means that the ground state orientation for a pair of he- 
lices can change drastically as flexibility is increased. 

Although we did not consider the general problem of 
discrete frustrations in aggregate systems, the two appli- 
cations wc covered demonstrate that such studies may 
be carried out easily on a case by case basis. For 13/6 
charges per turn F-actin paracrystals, our findings ap- 
pear to indicate that the observed twistings in aggre- 
gate cannot be explained by the electrostatics of the pair 
interaction. Given the recently presented results indi- 
cating that twisting may limit bundle width in protein 
linked F-actin aggregates [11], it seems plausible that lo- 
cal twisting may play some key role in the mechanics of 
counterion-induccd aggregation of F-actin as well. For 
A-DNA, our brief consideration of a hexagonally packed 
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system indicated that the dominating rational term com- 
pletely averages out for this geometry. This further 
strengthens the plausibility that non-hexagonal packing 
structures are observed for A-DNA to reduce this frus- 
tration. 

We conclude with a few comments on local twisting. 
First, wc note that an irrational-rational transition is still 
expected for flexible systems, though the boundaries of 
the phase diagram may be changed dramatically depend- 
ing on the stiffness of the system. As shown above, an 
infinitesimal twisting to a rational system can allow for 
a finite energy gain while costing a negligible amount of 
torsional twisting energy. Thus rational transitions are 
still expected to often occur. In highly flexible systems, 
however, local discommensuration states are allowed and 
compete with the rational, fully in register, states. At fi- 
nite temperature, the entropic gain associated with such 
states may make them more favorable. Second, wc reit- 
erate that the pair interaction studied in this paper may 
not accurately model highly flexible systems. For ag- 
gregates of DNA, we expect the pair interaction to be 
accurate, as DNA is known to have an unusually high 
torsional rigidity [l^l . Actin on the other hand is known 
to have an especially low torsional rigidity , implying 
that corrections to the pair interaction might be signifi- 
cant in this C&rSS, ClS mentioned above. 
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APPENDIX: THE COULOMB SUM 

In this appendix we outline the steps taken to express 
the perfect helix interaction energy in a form which may 
be easily evaluated numerically. Formally, the energy is 
given by Eq. ([T]), which we rewrite below. 



E 



E 



exp[— flgi?] 
R ■ 



(A.l) 



The sum here is over all pairs of charges, one taken from 
each helix, where, from Eq. R is given by 



R^ = cf - 2rd[cos{ei -cj))- cos(6'2 - </>)] 



-2r^[l - cos(e'i - 



-O2 + C\ 



(A.2) 



To begin we apply the identity 

r'^/'^ exp[iA:as - (fc^ -I- R^)t]dt, 



exp[-asi?] _i,2 

= TT ' 



R 



and then replace the exponentiated cosines in (jA.ip using 
the series representation 



exp[ycos6'] ^ Ii{y) cyi^[iie]. 



(A.4) 



Here /; is a modified Bessel function of the first kind. 
Then, the expression that will yield the interaction en- 
ergy is 

TT-^'^ / t-^'^Ii,{2drt)Ii,{~2drt)Ii,{2rh) 



exp 



At 



t[d^ + 2r^ 



-c? 



+ih{di -(f)) + ii2{e2 - 0) + ^^3(^1 - 02) 



dt. 



(A.5) 



We now let ni and 712 in Eqs. ^ and ^ be given by 



ni = TOa + mb 
n2 = mb. 



(A.6) 
(A.7) 



The sum on m^, may then be separated out. It is 
exp[^(/i + l2)mtip] 

TTif, — — 00 

00 

= 2n S{^lj{li+l2)~2TTk). (A.8) 

fc— — 00 

If 'ip/2n is irrational; then the only possible way in which 
an argument of one of the delta functions on the right 
hand side of Eq. (|A.8[) can be zero is if h + h = 0. 
For the time being, we will assume that this is the case. 
Then, we have an infinite contribution from all terms of 
the form li = —l2- To see what happens in this case, 
we set li = —I2 on the left hand side of Eq. (jA.8|) . The 
result is that we are summing over one for each charge 
on one of the helices. The reason that we end up with an 
infinite result is the sum as defined in Eq. (jA.8[) contains 
an infinite number of terms. If we are interested in the 
energy per charge, we simply take one of them, having 
set li + I2 = 0. This leaves us with the following sets of 
sums and integrals to perform 



E 



dt {^t)-^/^Ii, {2rdt)I_i, i-2rdt)Ii, {2rH) 



exp 



4i 



t{d^ + 2r^) - ahiruai' - At/. + 



+i{h + l3){ma'4' - AiP) 



(A.9) 



We now split the range of integration over t into two 
parts. When t > l/{aip)'^ there is reasonably rapid damp- 
(A.3) ing of the sum over ma- Let's call the lower limit for this 
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range of t values T. For the upper range of t values, we 
invert the sum over to get 

oc 

J2 /i3(2r2t)e*'^(™»^'^'^) = exp[2r2tcos(ma - AV')]. 

l3 = -00 

(A.IO) 

Through application of the identity 



Inix) = 



2tt 



^—in9-\-x cos £ 



(A.ll) 



the sum on li may be evaluated to give 

oc 
/i— — oo 



1 r'' 

— / exp[2r(ii(cos 6* 
27^ 7o 

- cos(0 - maV" + AV- + 2T:k)]de 



(A.12) 



Reconstructing the expression, we are left with the fol- 
lowing integral and double sum 



1 °° 

— y 



27r 



m a — — oo 
J2 , 



de I dt (7rt)-^/2exp 



t{d'' + 2r^)- a^iphirua - [AiP - Q/ij^f + 2rdt[cos6' 
— cos(0 — TOaV' + ^V')] + 2r'^f cos{inaip — Alp) . 

(A.13) 

It is important to note that because of the damping of 
the sum over TOq the number of terms actually summed 
over is not large. Furthermore, the integration over 9 is 
limited to a finite range. Finally, the integration over t 
is exponentially convergent. 

In the small-t regime, we concentrate on the sum over 
TOa in Eq. (jA.9[) . Here, we use the following identity, 
based on the Poisson sum formula 

CXD OO 

J2 f{m)= / /(m) ^ e^-'''''"'dm. (A.14) 

„ J — oo 



m— — OO 



fc' — — c 



This leads to the following expression to be summed and 
integrated over in the low-i regime. 

V / dt\li,{2drt)I^i,{-2drt)Ii,{2rH) 



exp 



^ - t{d^ + 2r^) 



(A.15) 



This last transformation has ensured that the sum on k' 
is now quickly damped in Eq. (|A.15[) . 



Things get a bit simpler when the interaction is un- 
screened, in that at least one of the integrals above can 
be expressed in terms of special functions. However, in 
the absence of screening the energy per charge diverges 
as expressed in Eq. (jA.ip . This divergence may be elim- 
inated if we subtract out the energy per charge for two 
interacting lines of charge. We may then add this term 
back on using the logarithmic expression for the energy, 
obtained by shifting the origin of the first line's potential. 
The term we subtract out is 



— / i ^exp 
aV Jo 



-d^ t 



dt. 



(A.16) 



which diverges in the small-t limit. The appropriate ar- 
rangement of canceling terms is the following: 



C(d,r,r) 



.2 ,,^2^ 



^1 -lh{2drtflo{2rH)c-K^[~t{d'' 
aip Jo t I 

- exp[-d^t]\ dt. (A.17) 



The simplification that results from removing the 
screening is in the integration over the upper limit. What 
we can do is express the integral over t in terms of an ex- 
ponential integral or of an error function. The relevant 
results are 



t 



t-^'^e-^'dt 



-dt = -Ei(-Tw) 



0Ferfc(VT^) 



(A.18) 



(A.19) 



The final expression for the interaction between two 
irrational helical arrays of charge, for the integration from 
to T is 



I.i,{~2drt)Ii,{2rh)e^p 

(2^fc' + (/i+/3W 



t(d2 + 2r^) 



Aa^ili'^t 

-C{d,r,T)-^\ogd. 
aip 



(A.20) 
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And for the integration from T to oo, we have 



where 



oo 

- F 



ma — — oo 



-2rd[cos 9 — cos(6' — niatp + Ai/;)] 



-2r^ cos(too'0 — 



erfc('ri/2y2 ^2r2 



)2 - 2rd[cos ( 



2Trl^ 



(A.24) 



- cos 



{6 - maip + AiP)] - 2r^ cos{maip - A^|:)]^/'^)d0 



Ei{-d^T) 

aip 



(A.21) 



Although the expressions in Eq. (jA.2ip are a bit cumber- 
some, their numerical evaluation, at least in Mathemat- Finally, the new version of Eq. (|A.15p is 
ica, is considerably easier. 

For rational tb = 2tt^^, the delta function generated 

by the summation over mt in Eq. (jA.SP is a bit more 
general. Now the condition is 



li + I2 = 2'iik^n^ 



(A.22) 



where is also an integer. After a bit of work, we find 
for the new version of Eq. (|A.2ip 



00 r 



+ 2r2 + a2^2(^^ _ _ ^)/^)2 



00 00 00 00 

E E E E 

/i — — 00 /3 — — 00 k— — oo k' — — oo 

r „2 

/,3(2r2t) exp 



dt— (2drt)/_,,+fc„^ {-2rdt) 



-^^-^(d^ + 2.^)-iM:±^l±W 



-2rd[cos 9(TOa, l^,(t>) - cos{Q {ma, 1 4,, 4>) ~ rua-ip + Alp)] 

-1-1/2 

— 2r^ cos(maV' ~ A?/') 

erfc(Ti/2[d2 + + a^ij^{ma - (At/; - O/i/jf 
-2rd[cosQ{ma,l^,4)) - cos(8(ma, ^0, </)) - maip + Aip)] 



-Mil + h) + i— -2Trk' + ikn^iAip - 



(A.25) 



-2r^ cos(matp - AtP)]^/^) + — F.i{-d^T), 



(A.23) 
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